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Abstract 

Self-gravitating systems have acquired growing interest in statistical mechanics, due 
to the peculiarities of the 1/r potential. Indeed, the usual approach of statistical 
mechanics cannot be applied to a system of many point particles interacting with the 
Newtonian potential, because of (i) the long range nature of the 1/r potential and of 
(ii) the divergence at the origin. We study numerically the evolutionary behavior of 
self-gravitating systems with periodical boundary conditions, starting from simple 
initial conditions. We do not consider in the simulations additional effects as the 
(cosmological) metric expansion and/or sophisticated initial conditions, since we 
are interested whether and how gravity by itself can produce clustered structures. 
We are able to identify well defined correlation properties during the evolution of 
the system, which seem to show a well defined thermodynamic limit, as opposed 
to the properties of the "equilibrium state" . Gravity-induced clustering also shows 
interesting self-similar characteristics. 
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1 Introduction 

Many-body systems with long range interactions are widespread in nature. 
In particular, gravitating systems dominate in an impressive range of spatial 
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scales (10 _1 joc— 10 8 pc) . It is very difficult to trace the specific role of gravity in 
the evolution and the statistical properties of such systems, since usually many 
different physical processes should be considered. For this reason, many theo- 
retical and numerical studies have focussed on the specific physical problem, 
rather than on the general statistical features of an infinite gravitating system. 
On the other hand, the usual statistical and thermodynamical approaches to 
the understanding of an infinite system with an interaction of some kind fail, 
or have to be taken very carefully, in the case of long range interactions. One 
of the properties of such interactions is that they make the system non exten- 
sive, that is, the system cannot be considered as made of smaller independent 
parts. Therefore one cannot study the typical properties of finite sample to 
extrapolate safely the properties of the whole infinite sample, as one does with 
short range systems. Many simple idealized examples have been considered to 
analyse the effect of non extensivity (e.g. [1-4]). They have given interesting 
clues on the behavior of non extensive system. On the other hand, it would be 
more satisfying to deal with three dimensional systems with a real interaction 
potential, as gravity, and a full treatment of the dynamics of the system. Our 
approach has been to deal with a system as realistic as possible, i.e. three 
dimensional N-body gravitating systems 



2 The simulations 

The model we investigate is usually a system with particles placed randomly 
in a cube with periodical boundary conditions. The particles have the same 
mass and are initially at rest. Since we want to understand the behaviour of 
an infinite system, we run a series of simulations with different number of 
particles, but with the same density. It is important to stress that, due to 
long range nature of gravitational interaction, one has to use subtle analytic 
tools to evaluate the force on a particle due to an infinite system. We use 
Ewald formula for this pourpose [5]. It is also important to observe that the 
gravitational potential in such systems is not well defined, since the potential 
energy of a particle would diverge. Neverthless, due to isotropy of the system 
on large scales, the force is well defined, since contributions from far away are 
balanced. We use a numerical code [6] which solves the iV coupled equations 
of motions in discrete time steps. A well known algorithm (tree code) is used 
to evaluate the forces in an efficient approximated way. The time integrator 
is a leapfrog (Verlet) algorithm, with individual time steps. For numerical 
reasons, the potential is smoothed at very small scales (much smaller than the 
mean interparticle distance), that is the divergence in the origin is removed 
and replaced with a smoother behavior. We have checked by changing the 
smoothing length that such procedure does not affect our results. 

We consider in the following the results of a series of simulations, with the 
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Fig. 2. Evolution of the kinetic energy per particle for simulations with 8 3 , 16 3 (filled 
circles), 32 3 (empty circles) particles. 

same number density but different number of particles, (8 3 , 16 3 , 32 3 ). In fig.l 
we show a picture of the evolution of particle positions in time. Note that 
the systems develop clusters on small scale first, then on larger and larger 
scales. When the typical size of the clusters is of the order of the box size, 
the evolution is influenced by finite size effect. Actually, even if the periodic 
system is infinite, the maximum number of particles in the simulation volume 
is fixed. 

In fig. 2 we plot the kinetic energy per particle versus time. This figure shows 
that the evolution of the kinetic energy per particle is independent of the 
number of particles in the system, in spite of the non extensivity, up to a time 
which increases approximately with InN . After this time, the kinetic energy 
becomes constant. This "equilibrium state", though, is a finite size effect. 
Therefore we argue that in the "thermodynamic limit" (defined as the limit 
behavior when N, V — * oo keeping the density constant) the behavior of the 
kinetic energy will be the same. It is also interesting to look at the distribution 
of velocities along an axis at a fixed time, shown in fig. 3. The picture refers to a 
given time before the finite size effects become relevant. A gaussian fit is shown 
for the reference. The distribution is non gaussian, but this is not surprising, 
since the system is not in equilibrium. It can be interesting to try to apply 
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Fig. 3. Normalized distribution of velocities along an axis for a 32 3 particle sim- 
ulation. The inner continuous line is the gaussian fit, the outer one is a stretched 
exponential fit. 

Tsallis' temptative velocity distribution (see e.g. [7]) on such data. We have 
also analysed clustering by means of the correlation functions. In particular we 
use the integrated conditional density T*(r, t) [8] and the integrated two point 
correlation function £(r, t) . Operatively they can be defined respectively as 
the mean average density and the mean average density fluctuation in a sphere 
of radius r centered on a particle of the system. Therefore a system is clustered 
on a scale ro if T*(ro,t) >> no, where no is the mean number density. The 
evolution of T*(r, t) is shown in fig.4 for simulations with different number of 
particles. Two important conclusions can be drawn from fig.4. First of all, the 
growth of correlations is approximately the same for all the simulations, in the 
regime where clustering takes places on scales small compared to the box size. 
Together with fig.2, this result suggests that the dynamics in this regime is the 
same for any number of particles, and therefore also in the "thermodynamic 
limit" . We can also reasonably infer that the infinite system has no equilibrium 
state. Another important remark is that the functional form of the T*(r,t) on 
the largest scales on which clustering is taking place is the same at all times. 
So we suppose that the clustering process could be explained in terms of the 
same dynamical processes, acting as time goes on on larger and larger scales. 



In two papers of ours [9,10] we have tried to explain the whole evolution of 
the correlations in such systems in terms of a discrete dynamics acting on 
larger and larger scales. At the beginning, the basic elements for the dynamics 
we consider are particles, which form clusters. At later times, the interactions 
between such clusters are mainly responsable for the further evolution of the 
system, and so on. This is a sort of renormalization of the dynamics. The 
justification behind such model lies on the fact that on scales much larger 
than the clustering scales the system is approximately isotropic, therefore the 
force due to far away objects almost balances. On the contrary, on the scale 
of clustering the system is anisotropic, therefore the force due to near clusters 
is large. In particular, it will be mostly in the direction of the nearest cluster. 
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Fig. 4. Evolution of T*(r, t) for simulations with 8 3 , 16 3 (lines with pluses), 32 3 (lines 
with filled circles) particles. Time flows from left to right. 
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Fig. 5. Comparison between £ evolution for a simulation with 32 3 particles (filled 
circles lines) and our theoretical prediction (smooth lines). The agreement for times 
smaller than the typical crossing time is quite good. 

This model is successful in providing: ( i ) an analytical expression for early 
time evolution of correlations, with no free parameter. The agreement with 
numerical experiments is given in fig. 5; ( ii ) a qualitative prediction for the 
rate of the clustering process at any time. 



3 Conclusions 

We analyse the problem of the evolution of an infinite system with long range 
(gravitational) interaction. We have studied it by numerical N-body simula- 
tions of periodical systems with different sizes. We identified some properties 
which can be inferred for the infinite system and we have given a model that 
could explain most of them. The relevant points of our approach are: 
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• We consider a thermodynamic approach to be highly problematic in the case 
of infinite self gravitating systems, since they presumably have no equilib- 
rium state and they are non-extensive 

• We propose a fully dynamical model instead which explains the out-of- 
equilibrium evolution of the system. The fundamental elements for the dy- 
namics in our approach are discrete objects. Discretisation turns out to be 
crucial to explain highly developed clustering, while standard (statistical 
mechanical or cosmological) approaches have mostly focussed on mean field 
or fluid like dynamics. In particular, the non-analiticity of the distribution 
of mass (that is, its discrete nature) is exported at larger and larger scale. 

• The evolution of clustering on a scale depends on what has happened on 
smaller scales. This is an opposite view with respect to the hierarchical 
clustering model, often considered in astrophysics. In particular, in such 
model clustering on a scale depends on initial conditions on that scale. 
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